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Abstract 

Iron oxides and oxyhydroxides are challenging to model computationally as competing phases 
may differ in formation energies by only several kJmol"^, they undergo magnetization transitions 
with temperature, their structures may contain partially occupied sites or long-range ordering of va- 
cancies, and some loose structures require proper description of weak interactions such as hydrogen 
bonding and dispersive forces. If structures and transformations are to be reliably predicted un- 
der different chemical conditions, each of these challenges must be overcome simultaneously, while 
preserving a high level of numerical accuracy and physical sophistication. Here we present compara- 
tive studies of structure, magnetization, and elasticity properties of iron oxides and oxyhydroxides 
using density functional theory calculations with plane-wave and locally-confined-atomic-orbital 
basis sets, which are implemented in VASP and SIESTA packages, respectively. We have selected 
hematite (a-Fe203), maghemite (7-Fe203), goethite (a-FeOOH), lepidocrocite (7-FeOOH), and 
magnetite (Fe304) as model systems from a total of 13 known iron oxides and oxyhydroxides; and 
use same convergence criteria and almost equivalent settings in order to make consistent compar- 
isons. Our results show both basis sets can reproduce the energetic stability and magnetic ordering, 
and are in agreement with experimental observations. There are advantages to choosing one basis 
set over the other, depending on the intended focus. In our case, we find the method using PW 
basis set most appropriate, and combine our results to construct the first phase diagram of iron 
oxides and oxyhydroxides in the space of competing chemical potentials, generated entirely from 
first principles. 

keywords: iron oxides and oxyhydroxides; phase diagram; density functional theory; modeling 
and simulation. 



I. INTRODUCTION 



Iron oxides and oxyhydroxides are abundant in nature; they are widespread in soils, waters 
and rocks, and are also found in living organisms, air dusts, meteorites, and Martian soils. 
p. 1-7] Iron oxides and oxyhydroxides have been the focus of numerous studies in the fields 
of geology, materials, soil, biology, and environmental sciences, and have broad applications 
in pigments, magnetic recording devices, medical imaging contrast agents, heavy metal se- 
questration absorbents. [l[ p. 2, 509-523] To date, 13 natural and synthetic iron oxides and 
oxyhydroxides (in addition to 2 hydroxides, see Ref. [l|, p. 2]) have been identified. These 
polymorphs have complicated structures (poor crystallization, ordering of vacancies, partial 
site occupancy), undergo a range of phase transformations, have characteristic magnetiza- 
tion states, and participate in a number of different types of interactions with contaminants 
and adsorbates. Given their ubiquity, it is surprising to find that the structures of some iron 
oxides and oxyhydroxides remain poorly understood, even after years of studies and nu- 
merous debates. In addition to this, size effects introduce a further complication, especially 
when we approach nanometer regimes, as shown in a recent review of the structure complex- 
ity, {2! in which the authors showed particle size, hydrous and hydrated environments, and 
synthesis processes all affect the observed structure. Collectively, these complicated issues 
have fueled constant interests in iron oxides and oxyhydroxides over the past decades. 

Like many materials, the development of characterization technologies and new samples 
often sparked renewed debates and led to new questions. One example is the debate on 
the origins of magnetite found in meteorites and magnetotactic bacteria. The magnetite 
nanocrystals from the Martian meteorite ALH84001 share many features with that from 
magnetosomes in terrestrial magnetotactic bacteria. [3] The similarities include unusual 
morphology, chemical purity, and crystallographic perfection. The similarities led to the 
proposal that the magnetite nanocrystals from the Martian meteorite were produced by 
biogenic processes, therefore provided strong evidence of lives in early Mars. 0, 5| This pro- 
posal was later questioned js, 7| and even dismissed js] because inorganic processes can also 
produce similar morphologies. However, the debate triggered new studies seeking reliable 
methods to identify origins or magnetite nanocrystals, and crystal size distributions 9| and 



oxygen isotope fractionation 
organic origins. 



10| have now been proposed to discriminate inorganic from 
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In recent years, computational modeling has opened up another potential way to solve the 
pending questions about iron oxides and oxyhydroxides. It generally requires electronic-level 
modeling methods to capture the rnagnetization states of iron oxides and oxyhydroxides, and 
density functional theory (DFT) ll| is able to solve electronic structures with desired ac- 
curacy at affordable computational cost. While DFT implementations have been routinely 
used to solve a wide range of problems in materials science, iron oxides and oxyhydroxides 
are particularly challenging for a number of reasons. Firstly, the energy differences between 
different solid phases or magnetization states may be as low as several kJmol^^, which is 
close to the resolutions of most DFT calculations, and necessitates energetic convergence cri- 
teria on the order of a few meV. Secondly, the underestimation of band gaps by DFT makes 
it difficult to depict the correct electronic structures of iron oxides and oxyhydroxides, most 
of which are semiconductors. A remedy to this problem is to include on-site Coulomb in- 
teraction to describe the strongly-correlated 3d electrons. [l2| Thirdly, the structures of 
iron oxides and oxyhydroxides may have partially occupied sites, or long-range ordering of 
vacancies (as in maghemite), which need large super cells and, accordingly, heavy compu- 
tation loads. Fourthly, the charge orderin g an d associated symmetry change in magnetite 



below the Verwey transition temperature 



13 



- |l9| are computationally intractable. Work- 



ing models proposed for charge ordering generally go beyond most DFT implementations. 
Fifthly, some iron oxyhydroxides have loose structure, eg. lepidocrocite (7-FeOOH), where 
the binding between layers relies on week hydrogen bonds and dispersive forces which, how- 
ever, are poorly described in DFT. And finally, various magnetization states in iron oxides 
and oxyhydroxides usually lead to slow convergence in calculations. 

The challenges of iron oxides and oxyhydroxides make computational modeling and sim- 
ulation non-trivial tasks, and work in this area tends to be sparse and sporadic. Despite the 



difficulties, DFT ca^ 
the past. [l3,ll8|,l2fl 



cu 



ations have been applied to some iron oxides and oxyhydroxides in 
31| These calculations incorporated different approximations, basis sets, 
and computational settings. It is therefore difficult to compare their accuracy and assess 
the methodology and algorithms, even though such a comparison is highly desirable for se- 
lecting computation tools in studying this difficult system. It also means that a systematic 
comparison between studies is not necessarily reliable, and cross-comparisons of different 
materials (such as those provided in phase diagrams) is not possible. However, when we 
seek to overcome this problem, we are confronted with the question of which is the most 
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appropriate technique to employ. 

In this study, we will present a comparative study between two implementations of DFT 
in calculating thermodynamic, magnetic and elastic properties of iron oxides, and assess the 
efficiency, accuracy and convergence, based on five iron oxides and oxyhydroxides, including 
hematite (Q;-Fe203), maghemite (7-Fe203), goethite (a-FeOOH), lepidocrocite (7-FeOOH), 
and magnetite (Fe304). Based on this large and consistent set of results, we are in a 
position to present the first environmentally sensitive phase diagram of iron oxides and 
oxyhydroxides, generated entirely from first principles, for predicting the thermodynamically 
stable structure as a function of the supersaturation of oxygen and/or hydrogen. 



II. COMPUTATIONAL METHODS 



A major difference among implementations of DFT is the choice of basis sets to expand 
the state space. Electronic wave functions can be constructed by linear combination of 
delocalized plane waves (PW's), or locally-confined atomic orbitals (LCAO's). The two 



basis sets have their advantages and disadvantages. [32|, |33| PW's have definite mathematical 
forms, are easy to implement, and have systematic convergence over cutoff energies, but 
their delocalized nature prevents linear scaling with the system size. LCAO's are flexible in 
terms of shape, size, and range, require much less number of orbitals compared to PW's, 
are localized thus suitable for spatial partition and linear scaling algorithms, but lack a 
systematic convergence and require extra effort to tune the LCAO parameters. We choose 



the DFT implementations in VASP (Vienna ab initial simulation package) 3J, |35| for the 



PW basis set, and SIESTA (Spanish initiative for electronic simulations with thousands of 



atoms) [361, |37| for the LCAO basis set. 



It is a known failure of local density approximation (LDA) or local spin density ap- 
proximation (LSDA) to accurately predict the ground state of bulk iron, while generalized 
gradient approximation (GGA) can reproduce the ferromagnetic BCC (body-centered cubic) 



ground state, 
and Ernzerhof 



Therefore, in our study, we choose GGA (in the form of Perdew, Burke 
39| ) for describing electron-electron interactions. For consistency, we use 



the same exchange-correlation functionals for both PW- and LCAO-based implementations, 
thereby enabling a complementary and detailed comparison between the two basis sets to 
assist others in this fleld. 
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A. Pseudopotentials 



In this study, we use pseudopotentials to describe core electrons and nuclei. For the PW 
basis set, we use the projector augmented wave (PAW) potentials from the pseudopotential 
libraries shipped with VASP. The reference states of valence electrons for generating the 
pseudopotentials of Fe is ScfAs^. The core radii are 2.30 Bohr, 1.1 Bohr, and 1.52 Bohr for 
Fe, H, and O, respectively. Nonlinear core corrections are included for Fe with radius of 2.0 
Bohr. 

For Fe, the 3d electron orbital overlaps with 3s and 3p core orbitals in real space, and has 
a small core radius of approximately 0.7 Bohr. 4s and Ap orbitals extrude further away from 
the nucleus, with radii of approximately 2 Bohr. These different core radii make it difficult 
to assign a common cutoff to all the orbitals, due to the short core radius of 3d orbital (~0.7 
Bohr), which requires very large cutoff of plane- waves (about 11000 eV) to converge the 
energy in 3 meV/atom. [40| One practice to eliminate the difference in core radii is to include 
3s and 3p as semi-core states, in place of the 4s and 4p states respectively. In this way, the 
reference state is not neutral (one 4s electron or two 4s electrons are excluded, assuming the 
ground state is 3(i^4s^ or 3(i^4s^), which is acceptable under the pseudopotential scheme. It 
is therefore possible to generate high-quality pseudopotentials with small core radii of around 
0.6-0.9 Bohr. The hard pseudopotentials can accurately reproduce all-electron calculations 
to excited states, but are computationally demanding. However, it has been previously 
shown the gain in quality of the calculation is not apparent when semi-cores are included in 



Ti and Cu. 



41| In order to reduce the computation cost, settings of ~2 Bohr radii have been 



found to be good compromise between efficiency and cost. The soft pseudopotentials often 
produce acceptable results in calculating lattice parameters, magnetization and electronic 
structures. 

For the LCAO basis set, we generate norm-conserving pseudopotentials according to the 
revised scheme of TrouUier and Martins. 42| A potential generated with the reference valence 
state of 3(i^4s^core radii of 2.0 Bohr, and partial core radius of 0.7 Bohr was used in previous 



studies. 



43 



44j | Since we wish to compare with our PW calculations, we have chosen the 



same reference states for Fe {3d'^4s^). According to our convergence tests, the core radii are 
2.0 Bohr for Fe, 1.1 Bohr for O, and 0.8 Bohr for H, smaller than those core radii of the PAW 
potentials for the PW basis set. Nonlinear core corrections are included for both Fe and O. 
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We found the nonlinear core radius of 0.7 Bohr provides the best match between pseudocore 
electron density and all-electron core-electron density. The pseudopotential of Fe generated 
with the same configurations has been used in studies to structure and magnetic properties 



of iron. 



44 



45| The nonlinear core radius of O is 0.7 Bohr, which is same as that in ref. 46 |. 
We test the transferability of the pseudopotentials by comparing atomic energies of excited 
states from pseudopotentials and from all-electron calculations. 

It is important to point out that a different nonlinear core radius is used for Fe, and 
nonlinear core correction is excluded for O in PAW potentials. These differences reflects 
to the degrees of compromise between efficiency and transferability. Fortunately, the pro- 
vision of a standard pseudopotential database (by VASP) allows for considerable testing 
in a large variety of situations, and the norm-conserving pseudopotentials for LCAO have 
been tested in the above-mentioned references. Therefore, we are confident that both sets 
of pseudopotentials should represent core electrons of Fe, H, and O, and are adequate for 
this comparative study. 

B. Basis sets 

PW's have a definite mathematical formula with no adjustable parameters. LCAO-based 
basis sets use the so-called pseudo-atomic orbitals (PAO's) whose shape, size, and range 
are configurable. The PAO's are mathematical functions with adjustable parameters, which 
must be optimized for specific systems, and the quality of the PAO's are critical to the 
simulation results of LCAO basis sets. In the present study we have optimized our PAO's 
by comparing simulated and known properties of simple structures, specifically, the lattice 
parameter of bulk body-centered cubic (BCC) Fe, and the bond lengths of H2 and O2 
molecules. 

The PAO's in the present study (for all the three elements) are of double-C plus polar- 
ization (DZP). The dimensionless parameter split-norm, which determines the splitting of 
different C functions, was set to 0.28 for Fe, 0.24 for O, and 0.65 for H. The large split-norm 
of H is in accordance with the large variation in the effective spatial extent of hydrogen in 
charged states. A similar value of 0.5 was employed during a study of the pressure effects on 
hydrogen bonds as reported in Ref. |47|]. Soft confinement has been applied according to the 



scheme proposed in Ref. 48| to avoid discontinuity of the functions at the cutoff distance. 
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TABLE I: PAO parameters of Fe, O and H. The and rc2are radii of double-C and 
polarization orbitals; V is soft-confinement potential, and is inner radius of the soft 
confinement; "Polar." represents polarization orbitals. 





rci(Bohr) 


(Bohr) V (Ry) n 


(Bohr) 


Fe 3d 


4.229 


2.292 


50 


3.81 


Fe 4s 


6.800 


5.363 


150 


6.12 


Fe Polar. 


6.800 




150 


6.12 


2s 


5.000 


2.580 







2p 


6.500 


2.497 







polar. 


3.923 




104.3 


0.00 


H Is 


4.971 


1.771 


2.07 


0.00 


H polar. 


4.988 




0.89 


0.00 



The parameters for generating the PAO's are summarized in Table [H In the DZP scheme, 
the numbers of PAO's per atom are 17 for Fe, 13 for O, and 3 for H. The results of bulk Fe 
and the gas molecules (H2, O2, and H2O) used to construct our basis sets, are provided in 
the Section UTTl 



C. GGA+C/ parameter 



The strong correlation effects of iron 3d electrons lead to splitting of d bands. Depending 
on the relative positions of oxygen 2p and iron 3d orbitals in valence bands, iron oxides and 
oxyhydroxides may be semiconducting or metallic, [l, p. 115-117] Both GGA and LDA tend 
to over-delocalize electrons and underestimate correlation effects and band gaps. Model 



Hamiltonian approaches are often used in such strongly correlated systems. [49j In these 
models, electrons hopping between atoms experience the effective Coulomb interaction U, 
which is defined as the energy cost for moving an electron between two atoms that both 
initially had the same number of electrons, or U = + — 2i?„, where En is the energy 
of an atom with n 3d (for transition metals) or 4/ (for rare earth elements) electrons. 49 1 
This energy fluctuations result in the formation of band gaps. As implementation of the 
model Hamiltonian approaches in DFT, the LDA+J7 (or GGA+U) method 12|, |50| includes 
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on-site Coulomb interactions among strongly correlated electrons. 

We point out that there exist alternative approaches to solve or alleviate the band- 
gap problem o 



correction 
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DFT, including hybrid HF-DFT functionals |5ll . |52|, and self- interaction 



54j |. These approaches (including the aforementioned DFT+f/) are being 



extensively tested in a large variety of chemical environments and becoming widely imple- 
mented. Particularly for strongly correlated systems, hybrid functionals have been shown 



to properly describe the magnetic coupling in and band gaps of NiO |55|, U O2 56| . Ce02 
and Ce203 [s^], plutonium oxides 58 1, and several strongly correlated solids 59 1. Various 
hybrid functionals have been developed and actively tested with other hybrid functionals 
and pure functionals. 



60 



range-separated hybrids [59| 



Among the recent developments of hybrid functionals, 
62|-|64| and Heyd-Scuseria-Ernzerhof hybrid functional 65 



the 



68| 



are very promising in tackling the correlation effects in solids. Hybrid functionals often give 
acceptable thermochemical results owing partly to their semi-empirical nature and the fit- 
ting procedure (such that the amount of exact exchange can be tuned to fit known physical 
and chemical properties). With the increasingly available options, it is, however, desirable to 
select those density functional approximations of nonempirical constraint satisfactory with 



least fitting parameters. 



69| Choices of the approaches may depend on the availability of 



implementations or computational cost. In this study, we have chosen DFT+U to account 
for the band-gap problems of DFT because it is implemented in both computation packages 
(VASP and SIESTA) which are analyzing herein. There are, of course, many other com- 
putation packages that use PW and LCAO basis sets with different compromise between 
accuracy and computation cost. Both the two packages we use for this study have users of 
broad interests, ranging from physics, chemistry, materials science, and biology. They thus 
should serve as robust computational tools for our study. 

For the strongly-correlated systems, the improvements to band-structure calculations pro- 
vided by DFT+U are substantial. 50| To demonstrate this, we test the GGA+U methods in 
the calculations of bulk Fe, iron oxides and oxyhydroxides in both PW and LCAO methods. 
Both packages have implemented the GGA+U method based on a simplified rotationally 
invariant formulation by Dudarev et al. [50]. In this implementation, only the effective 
Coulomb repulsion U^s = U — J is significant. In our study, the on-site Coulomb interac- 
tions are included only for strongly correlated Fe 3d electrons, but not for the electrons of 
O or H, or other types of electrons of Fe. 



8 



The DFT+[/ method has previously been employed to study magnetite 



hematite l2 



22 



17 



18|, 



70|, goethite |29|, and maghemite 31|, for which the parameter U varies 
between 2 eV to 5 eV. Cococcioni and Gironcoli suggested U^s of bulk iron to be ~ 2.2 eV us- 
ing a linear-response approach; [71] Anisimov and Gunarsson gave rather large Ucs of about 



6eV. 



49| Rollmann et al recommended Ues of 3.0 eV through their study of the electronic 



structure of hematite. 



Punkkinen et al suggested a much smaller value (~ 1.0 eV) for 
hematite, designed to reproduce experimentally observed features of the electronic structure, 
such as the crystal field induced band splitting, 
in the study of vacancy ordering of maghemite. 



20| Grau-Crespo et al. used Ucs = 4.0 eV 
jsi ] The difference may originate from im- 
plementations of the DFT+^7 method, pseudo-core configurations, and even basis sets. In 
the present study, the parameter of Ues is chosen so that the calculated band gaps and 
lattice parameters both match the experimental values. We found Ues = 4.5 eV provides 
best match the experimental band gaps of hematite (see section IIII Cp and goethite (not 
shown), and acceptable lattice parameters of all the iron oxides. The same Ucs was used for 
all the iron oxides for consistency. 

It is worth noting that standard DFT (LDA or GGA) reproduces thermodynamic prop- 
erties very well, sometimes exceeding the predictions of the DFT+f/ method in comparison 
to experiments. However, DFT+[/ methods provide much more accurate predictions of 
electronic structures. Ideally, first principle methods should accurately predict both ther- 
modynamic and electronic properties, but this remains a goal of those involved in the devel- 
opment of new density functionals. In our calculations, we compare the results of GGA+U 
with GGA, to assist others in selecting the most appropriate approach for their work. 



D. Computational Settings 

To facilitate a cross-comparison, we have used consistent settings for all the iron oxides 
in both basis sets. The fe-points for sampling over the Brillouin zone were generated using 



Monkhorst-Pack scheme. [7^ For a primitive cell of BCC Fe, a A;-mesh grid of 23 x 23 x 23, 
which corresponds to 364 irreducible fc-points in the first Brillouin zone, can achieve con- 
vergence of total-energy in 2 meV/atom when using the PW basis set. For the LCAO basis 
set, a A;-grid ofl7xl7xl7 can reach the same convergence of energy, and the number of 
A;-points is 2457. One immediately notices the large difference in the numbers of fc-points 
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TABLE II: Sizes of /c-meshes and numbers of fc-points in the calculations using the PW 
and LCAO basis sets. The numbers in the column of PW are the numbers of irreducible 
A;-points in the first Brillouin zone, and the numbers in the column of LCAO are the 

numbers of trimmed fc-points. 



/c-grid PW LCAO 



magnetite 


4x4x4 


10 


44 


hematite 


4x4x4 


13 


64 


maghemite 


2x2x2 


1 


8 


goethite 


4x6x4 


24 


60 


lepidocrocite 


8x4x8 


32 


150 



in the PW and LCAO basis sets. This is due to the different symmetrization treatments 
in the two programs. The VASP code utilizes crystal symmetries to calculate the charge 
density, forces, and stresses. The symmetry elements of the crystal structure greatly reduce 
the number of necessary A;-points for adequate sampling. The SIESTA code is designed 
for large systems, as its name indicates, and symmetry constraints are usually excluded. 
SIESTA only trims a small amount of redundant /c-points from the constructed grid. Alter- 
natively, SIESTA uses molecular dynamics (MD) algorithms for geometry optimization over 
an auxiliary supercell. This difference in symmetrization leads to a very different number 
of /c-points used in sampling the band energies, however, the convergence criteria of fc-mesh 
density with respect to total-energies are set to 1-2 meV/atom in both basis sets. The sizes 
of fc-mesh and numbers of fc-points used in the calculations are shown in Table HIl 

In the PW basis set, we find a plane- wave cutoff of 800 eV can achieve convergence in 
the total-energies to below 1.0 meV/atom for all the five iron oxides and oxyhydroxides 
considered in our study. For bulk Fe, a smaller cutoff (600 eV) is able to achieve the 
same convergence. In the calculations of isolated O2, H2, and II2O molecules, the PW 
cutoffs are 850, 600, and 850 eV, respectively. With these cutoffs, the difference in total- 
energies can be reduced to less than 1 meV/atom, which is the limiting resolution of the 
DFT implementation. SIESTA uses a finite real-space grid over which integrations are 
performed to calculated energies, forces, stresses, and dipoles. The fineness of this finite 
grid is determined by a cutoff value, which is equivalent to the PW cutoff in the PW basis 
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set. There are subtle differences between these equivalent settings across the two basis sets. 
In the LCAO basis set, wave functions are constructed using atomic orbitals, and the cutoff 
should only affect the accuracy of integration; in the PW basis sets, the plane waves are 
also used to construct the valence wave functions, so the cutoff has a larger impact on the 
quality of calculations. After the convergence tests, we chose cutoffs of 5130 eV for bulk Fe, 
4080 eV for O2, 2040 eV for H2, and 6800 eV for all the 5 iron oxides and oxyhydroxides, so 
that the total-energies converge below 2 meV/atom. 

Gaussian (in the PW basis set) or Fermi-Dirac (in the LCAO basis set) distribution 
functions are used for electronic occupations for the molecules (H2, O2, and H2O), iron 
oxides and oxyhydroxides; Mehfessel-Paxton functions of order 1 is used for bulk iron. The 
smearing width or electronic temperature has been set to 0.03 eV for all the iron oxides and 
oxyhydroxides, and 0.05 eV for bulk Fe, in both basis sets; relatively small value (0.15 eV) 
is used for the isolated O2 molecule, and large values (0.4 or 0.5 eV) are used for the isolated 
H2 and H2O molecules (in both basis sets). These values are chosen so that the energies 
diverge by less than 2 meV/atom compared with smaller smearing widths. 

Energy minimizations to all the structures are conducted using conjugate gradient (CG) 
algorithms with the force convergence of 0.005 eV/A. For the iron oxides and oxyhydroxides, 
geometry optimizations of unit cells are done with a convergence criterion of 0.005 GPa for 
the stress tensor components. For the isolated molecules (O2, H2 and H2O), a large super 
cell of 10 X 10 X 10 is used (in both basis sets). 

E. Magnetization states 

Iron oxides and oxyhydroxides undergo magnetic phase transitions at different tempera- 
tures. Most of them are antiferromagnetic or ferrimagnetic at temperatures below their Neel 
or Curie temperatures. Magnetite and maghemite are ferrimagnetic; hematite, goethite, and 
lepidocrocite are antiferromagnetic. P, p. 123] In this study, we consider alternative magne- 
tization states in addition to those observed experimentally. By comparing the energetic 
stability of different magnetization states, we are able to test the validity of our calcula- 
tions. In general, a non-spin polarized state, a ferromagnetic state, and several other initial 
spin-polarization states are included. However, we only consider collinear magnetization 
states, which are most commonly observed in iron oxides and oxyhydroxides at low temper- 
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atures. 

For consistency, we have also included spin polarizations when calculating the properties 
of the isolated molecules, even though H2 and H2O are non-magnetic (or diamagnetic) . 
In both PW and LCAO calculations, the net spin moments of H2 and H2O are zero, in 
agreement of experimental observations. The spin moment of O2 is 2.0 fiB using both PW 
and LCAO basis sets. 



F. Elastic properties 

In this study we calculated bulk moduli of each solid material by fitting to Birch- 
Murnagham equation of state. [73] In addition to this, we calculate the elasticity tensors of 
bulk Fe, iron oxides and oxyhydroxides using a finite-difference method. In this method, a 
series of strains are applied to the equilibrium unit cell, the total-energies of the strained 
structures are calculated, and the elasticity tensor components Cij are calculated through: 

where E is the total-energies of strained structures, Eq is the total-energy of equilibrium 
structure with zero stresses, and e is the applied strain. The subscripts i and j are of matrix 
notations, j?^ . p. 134] The strains are grouped into a number of transformations, which are 
chosen in accordance with the crystal symmetry of the structures. For each transformation, 
6 strains of ±0.015, ±0.010, and ±0.005, in addition to the equilibrium structure, are used 
for linear least-square fitting to calculate the tensor components. 

We developed a computer program to calculate elastic constants of crystals by using ab 
initio packages as backends. Since this method only requires total-energies, which can be 
calculated using many computation packages, we can make consistent comparisons by using 
the same strains. This method and program have been previously tested in calculating 



elastic constants of Co [75.] and Ni-B alloys [76|. 



III. RESULTS AND DISCUSSIONS 

In the following sections we will focus on presenting results of our detailed comparisons 
between the PW and LCAO basis sets, as well as the physical comparisons being made in 
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energetic stability, lattice parameters, elastic properties, and magnetization states of our 
collection of iron oxides and oxyhydroxides. 



A. Bulk Fe 

The ground state of bulk iron is of body-centered cubic (BCC) structure (space group 
ImSm, No. 229) and ferromagnetic. Fe is a well-behaved system within the framework of 
standard DFT-GGA, which predicts correctly thermodynamic properties such as energetic 
stability and lattice parameters. With on-site Coulomb interactions, first-principles calcu- 
lations can improve the predictions to electronic band structures. As mentioned above, the 



parameter Ues may 



interpretations. 



49 



vary 



71 



Tom below 1 eV to about 6 eV, depending on the methods and 



77H79l| As our focus is on thermodynamic properties, we apply mild 
on-site Coulomb interactions with Ues = 1.0 eV when calculating the properties of bulk 
iron. We choose this value of Ues because it improves the predictions of the lattice constant 
and cohesive energy in PW basis set (see Fig. [Hand Table [XTTll . In general, we find that the 
lattice constant of Fe increases almost linearly with Ueg. This is because on-site interactions 
alter charge density around Fe atoms, weakening the metallic bonding strength, similar to 
that observed in NiO. 50|] The spin moment, which is sensitive to changes of atomic volume, 
also increases with Ues- 

The calculation results are summarized in Table lllli Both calculations using PW and 
LCAO basis sets reproduce experimental lattice constants within 1.5%. All the calculations 
overestimate cohesive energy with respect to experimental measurement. LCAO overesti- 
mates by about 1.5 eV. The difference between calculation and experiment is much smaller 
in PW basis set. GGA using PW overestimates by about 0.5 eV, while GGA+f/ reduce 
the overestimation to about 0.1 eV. The calculations of the spin polarization moments (not 
including orbital moments) compare favorably with experiments, at around 2.5 /i^. 

For the calculations of the bulk moduli, GGA+U using PW best matches the results from 
experiments, while we find other methods overestimate the values by between 8% (GGA with 
LCAO) to 18% (GGA with PW). Both fitting errors and temperature effect may contribute 
to the difference between calculations and experiments, because bulk moduli are calculated 
at ground state from Birch-Murnagham equation of state, and experiments are conducted 
at the thermodynamic standard state. For the calculations of the elasticity tensors, we 
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FIG. 1: Dependence of (a) lattice parameter and (b) spin moment on Ucs in PW-based 
calculations. Experimental values are marked with horizontal lines. 



find that GGA using LGAO provides the best overall results, and other methods either 
underestimate C44 or overestimate Cn and C12 significantly. In particular, the components 
of the elasticity tensor calculated using GGA and the PW basis set can be considerably 
different to the experimental values, especially in the case of cn, but they are very close 
to those in recent calculations using exact muffin-tin orbitals and PBE functionals. [sol 
However, GGA+f/ tend to considerably underestimate C44 using both PW and LGAO. The 
differences between calculations and experiments may include defects in single crystal Fe 
being measured, extrapolation to ground state, anharmonic effects, and numerical error in 
the calculations. 

As shown in Table IIIIl GGA+^7 generally offers a small improvement over GGA in 
calculating the lattice constant and cohesive energy of bulk iron, at the expense of apparent 
underestimation of C44. 



B. Gas molecules 



As stated above, we have calculated the binding energies and bond lengths of H2 and 
O2, and bond angles of H2O (see Table ITVl) . Except for the binding energy of O2, the 
calculation results match experimental values within 2.5%. The significant overestimation 
of binding energy of oxygen dimer (and all other first row elements with more-than-half- 
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TABLE III: Calculation results of ferromagnetic Fe in BCC structure, oq: equilibrium 
lattice constant; Ec- cohesive energy; Bq: bulk modulus; M: spin polarization moment per 
Fe atom; Qj-: elasticity tensor components. Statistic errors of linear least square fitting are 



included for Cij. 



PW LCAO Expt. 



GGA GGA+C/ GGA GGA+C/ 



ao(A) 2.833 2.878 2.868 2.909 2.87(") 

Ec (eV) 4.94 4.36 5.97 5.64 4.28^^) 

M{pb) 2.20 2.79 2.31 2.67 (about 2.5) 

Bo(GPa) 198.4 164.6 182.1 188.0 168. 3^^) 

cii (GPa) 302.9±1.4 207.4±0.1 262.3±8.1 230.6±1.5 243.1 239.3^^), 297.8(-^) 
ci2 (GPa) 151.6±1.5 151.0±0.2 126.8±14.7 165. Oil. 8 138.1^^^), 135. 8^^), 141. 9(-^) 

C44 (GPa) 97.8±1.4 58.9±0.2 97.0±1.8 73.1±1.4 121. 9^'^), 120.7^^), 106. 7(-^) 

(a) Ref. [81, p.23]. 

(b) Ref. [82j. 

(c) Ref. \81, p.59]. 

(d) Ref. [83]. 

(e) Ref. [84]. 

(f) Ref. [80]. 

filled j9-orbitals) by DFT is due to the insufficient description of exchange energy, and lack 



of error canceling because of different electron shapes of O and O2. [85|, |86| Since the energy 
of oxygen dimer tend to be canceled out when calculating the energy differences between 
different phases, this overestimation is unproblematic in calculations of compounds, bulk 
iron oxides and oxy hydroxides. 



C. Hematite (a-Fe203) 

Hematite belongs to the trigonal space group of R3c (No. 167), and is isostructural with 
corundum AI2O3 or Ilmenite (FeTiOa). It is one of the most thermodynamically stable and 
abundant phases among all of the iron oxides and oxyhydroxides. ij, p. 6] Each rhombohedral 
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TABLE IV: Calculated properties of gas molecules. L stands for bond length, E for 
binding energy, for formation enthalpy, and a for bond angle of H-O-H in H2O. 



PW LCAO Expt 



(a) 



i^H-H() 


0.7500 0.7465 


0.7414 


E^_n (eV/bond) 


4.538 4.749 


4.521 


^o-o() 


1.2323 1.2422 


1.2074 


^0-0 (eV/bond) 


6.8074 6.2181 


5.1697 


^0-h() 


0.9575 0.9754 


0.9575 


aH-o-nO 


104.46 104.93 


104.51 


^f(H20) (kJmol-i) 


-243.8 -234.3 


-241.8 



(a) Ref. [87], pages 9-22, 9-24, 9-57, 9-58. 



unit cell contains 4 Fe atoms, distributed over 2 interlayer spaces of cation layers. Hematite 
is antiferromagnetic with all Fe ions in the same close-packing layer (perpendicular to the 
trigonal axis [0001]) having parallel spin moments, and different layers having antiparallel 
spin moments, noted as AFM (see Fig. [2j). At low temperatures below about 250 K, the spin 
moments change direction from perpendicular to parallel to the trigonal axis, keeping the 



antiferromagnetic configuration; 88| and no reports have found that the crystal structure 
changes at this magnetic transition. In order to validate our calculation results on hematite, 
we have included another two types of antiferromagnetic configurations in which Fe ions in 
the same layer have antiparallel spin (noted as AFM' and AFM" ; see Fig. [2]), a ferrimagnetic 
(FiM), a ferromagnetic (FoM), and a non-magnetic (NM) configurations. 

Ignoring on-site interactions leads to significant underestimation to the band gap (0.5 eV 
in calculation compared with 2.2 eV from experiment; see Fig. |3]). We find the calculated 
band gap linearly increase with ?7efj, as shown in Fig. |H In choosing the parameter U^s 
in GGA+[/ calculations, we fit the band gap to experimental value (about 2.2 eV [l]). To 
reproduce the experimental value of 2.2 eV, Ues should be between 4.0 eV to 5 eV. We have 
therefore adopted Ucs = 4.5 eV, and used this value consistently in our calculations to all 
the iron oxides and oxyhydroxides (in addition to hematite) using both PW and LCAO basis 
sets. 

The calculated thermodynamic and elastic properties for hematite are listed in Tables |V] 
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(a) AFM (b) AFM' (c) AFM" 

FIG. 2: Antiferromagnetic configurations of hematite. Antiparallel spin moments are in 
blue and cyan colors. Tlie rliomboliedral unit cell is marked by green lines (color online). 

and IVII We see both PW and LCAO correctly predict the lowest energy state of the 
antiferromagnetic configuration AFM, in agreement with experimental observations. The 
calculated lattice parameters match experimental values within 3% for this stable config- 
uration. In this case we find PW does a better job of reproducing the lattice parameters 
than LCAO. In both basis sets the lattice parameters from GGA+^7 are slightly larger than 
those from GGA. The largest difference between the two basis sets is the spin polarization 
moment of the metastable ferromagnetic state (FoM). In PW, the average spin moment 
is low (1.00 /ifi), in contrast with that in LCAO, where the average spin moment is high 
(3.42 /is). Both basis sets agree on the energetic order of the magnetization states, pre- 
dicting that AFM<FiM<FoM<NM (formation energy increasing), but differ about the two 
antiferromagnetic states (AFM' and AFM"). 

In the case of hematite, the calculated elastic constants using PW compare favorably 
with those determined using LCAO. Both basis sets give almost zero C14 (within numerical 
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FIG. 3: Density of states of antiferromagnetic (AFM) hematite. Up- and down-spin states 
are in solid and dashed hnes, respectively. The numbers are band gaps for f/eff = eV (no 
correction of on-site interaction) and Ues = 4.5 eV. The Fermi energy is shifted to eV. 
The density of states is calculated with fully relaxed structures. 



errors), while the calculation results of C44 differ by ~ 30 GPa between PW and LCAO. 
We notice that GGA+U produce larger values of C33 than GGA using both PW and LCAO 
basis sets, indicating on-site interactions strengthen bonding along the trigonal axis. All the 
elasticity tensors satisfy the elastic stability condition, which means hematite is elastically 
stable in all the four calculations. 



D. Maghemite (7-Fe203) 

Maghemite occurs as a weathering product of magnetite, and resembles magnetite in 
structure and magnetic properties. The Fe ions are all in the trivalent state, with balanc- 
ing vacancies to maintain charge neutrality. The crystal structure of maghemite has been 
characterized to be cubic, the same as magnetite, with partially occupied vacancies at oc- 



tahedral sites. 



89 



90| Depending on the ordering of cation vacancies, maghemite may be 



classified in either cubic (FdSm or P4332) or tetragonal (P4i2i2) space groups. Somogyvri et 
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FIG. 4: Dependence of band gap of hematite on Ues in GGA+U calculations. 



al reported long-range ordering of vacancies in powder neutron and XRD of nano crystalline 
needle-shaped maghemite, and classified maghemite to be in P4i2i2 space group. 90|] Using 



powder neutron diffraction, Greaves proposed the true symmetry of maghemite is tetrag- 
onal P432i2 instead of cubic P4332. SqI] The lattice parameters of the tetragonal cell are 
a = 8.3396 A, and c = 24.966 A which is slightly smaller 3a. 8^ Grau-Crespo et al sorted out 
the energetic order of various possible vaccancy ordering and found the tetragonal P4i2i2 
configuration has the much lower energy (by > 32kJmor^) than non-tetragonal configura- 
tions using classical interatomic potentials. 3l| The configurations P4i2i2 and P432i2 bear 
much similarity in structures and thus should have very similiar energetics. In this study 
we adopted the configuration proposed by Greaves [89| (tetragonal P432i2 symmetry), each 
unit cell having 160 (64 Fe and 96 O) atoms. Maghemite is ferrimagnetic below Curie tem- 
perature which is estimated to be between 820 K and 960 K. The Fe atoms at the tetrahedral 
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TABLE V: Calculated thermodynamic properties of hematite. The numbers in parenthesis 
are the errors in per cent compared with experimental values. The spin polarization is 
total spin moment of a unit cell (so that antiferromagnetization has exactly zero spin 
moment) divided by the number of Fe atoms. The energetically most stable magnetic 
states are in the first row of each basis set (PW or LCAO), GGA+[/ calculations start 
from the energetically most stable magnetic states. FoM stands for ferromagnetic, FiM for 
ferrimagnetic, AFM for antiferromagnetic, and NM for non-magnetic. GGA+U 
corresponds to the magnetization state of the lowest energy in GGA calculations. 





a (A) 


c(A) 


M(/iB) 


E{ (kJ mol"^) 


AFM 


5.005 (-0.6) 


13.884 


(+1.0) 


0.00 


-628.0 


AFM' 


4.841 (-3.8) 


13.183 


(-4.1) 


0.00 


-561.2 


AFM" 


5.044 (+0.2) 


13.850 


(+0.7) 


0.00 


-588.1 


PW FiM 


4.977 (-1.1) 


13.707 


(-0.3) 


1.50 


-592.3 


FoM 


4.783 (-5.0) 


13.333 


(-3.0) 


1.00 


-569.3 


NM 


4.733 (-6.0) 


13.511 


(-1.8) 




-543.0 


GGA+;7 5.074 (+0.8) 


13.874 


(+0.9) 


0.00 




AFM 


5.091 (+1.1) 


13.995 


(+1.8) 


0.00 


-676.8 


AFM' 


5.167 (+2.6) 


13.781 


(+0.2) 


0.00 


-640.1 


AFM" 


5.137 (+2.0) 


13.955 


(+1.4) 


0.00 


-638.5 


LCAO FiM 


5.056 (+0.4) 


13.809 


(+0.4) 


1.50 


-629.6 


FoM 


5.026 (-0.2) 


13.955 


(+1.5) 


3.42 


-585.6 


NM 


4.751 (-5.6) 


13.755 


(-0.0) 




-536.7 


GGA+U 5.183 (+3.0) 


14.072 


(+2.3) 


0.00 






5.034 


13.752 


0.00 


-823 ~ -828 



(a) Measured at room temperature and 0.1 MPa. Ref. [1, pp. 11 and 187] 



sites (where each Fe forms bonds with 4 nearest O atoms) have antiparallel spin moments 
with those at the octahedral sites (where each Fe forms bonds with 6 nearest O atoms). 
Specifically, the 40 Fe atoms in the 3x1 supercell at positions [1/8,5/8,0], [3/8,1/8,2/24], 
[1/8,7/8,2/24], [7/8,5/8,2/24], [3/8,3/8,0], [7/8,7/8,0] consist the majority spin component. 
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TABLE VI: Calculated bulk moduli and elastic constants of hematite. Unit: GPa. 





PW 


LCAO 




GGA GGA+U 


GGA 


GGA+U 


B 


174.4 190.3 


173.4 


176.1 


Cll 


325.0 ± 19.0 355.4 ± 13.9 


310.3 ± 12.1 319.6 ± 18.6 


Cl2 


131.8 ±5.6 132.1 ±6.4 


137.2 ±0.3 


125.6 ±4.3 


Cl3 


105.8 ±22.9 116.0 ±4.0 


114.6 ±6.9 


104.5 ± 1.5 


Ci4 


1.2 ±2.9 -5.6 ±6.4 


6.1 ±6.9 


5.5 ±9.4 


C33 


264.2 ±25.1 307.2 ±3.3 


255.6 ±3.8 


294.4 ±9.8 


C44 


103.0 ±6.4 110.6 ±7.0 


78.4 ±4.9 


80.1 ±5.2 



and the 24 Fe atoms at [4/8,6/8,1/24], [0,2/8,1/24], [2/8,4/8,3/24] consist the minority spin 
moment. Measurements of magnetic moments (spin polarization ± orbital moment) showed 
Fe atoms at the octahedral and tetrahedral sites have unequal spin moments: 3.54 jj,^ versus 
4.03 hb, [90] or 4.18 fiB versus 4.41 /xb. [89] In addition to this ferrimagnetization, we include 
ferromagnetic and non-magnetic states for validation of the calculation results. 

For the ferrimagnetic state, the calculated lattice constants match experimental values 
within 1.7% (see Table rvTTj) . with the exception of results from GGA±f/ using LCAO, in 
which the errors are about 3.1%. In general, the results from PW calculations are closer 
to those reported from experiments. If we omit spin polarization, the lattice parameters 
simultaneously decrease in both calculations, and the mismatch in lattice constants between 
the calculations and experiments increases to —4.4%. Including on-site interaction leads to 
a lattice expansion of about 1% in both PW and LCAO. Both calculations reproduce the 
correct magnetic ordering, predicting that that ferrimagnetic state has a lower formation 
energy than non-magnetic and ferromagnetic states. By comparing Tables IVIII and |Vl one 
immediately sees maghemite has a higher formation energy, and thus less thermodynamically 
stable than hematite. 

In the case of maghemite, the calculated elastic properties are very similar, both in trend 
and numbers, when we compare the PW and LCAO calculations (see Table IVIIII) . The 
diagonal components (cn and C33) of the elasticity tensor are noticeably larger in GGA+U 
than GGA; the shear moduli (C44 and C55) are also slightly larger when using GGA+U. We 
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TABLE VII: Calculated thermodynamic properties of maghemite. See Table I for 

annotations. 





a (A) c(A) M(/iB) 


El {kJ mol-i) 


FiM 


8.363 (+0.3) 25.034 (+0.3) 1.25 


-623.9 


PW FoM 


8.192 (-1.8) 24.563 (-1.6) 2.75 


-568.9 


NM 


8.026 (-3.8) 24.871 (-4.4) - 


-508.5 


GGA+U 


8.428 (+1.1) 25.237 (+1.1) 1.25 




FiM 


8.480 (+1.7) 25.374 (+1.6) 1.25 


-661.4 


LCAO FoM 






NM 


8.059 (-3.4) 24.153 (-3.3) - 


-496.5 


GGA+U 8.598 (+3.1) 25.718 (+3.1) 1.25 




Expt.(") 


8.34 24.97 


-806--813 



(a) Measured at room temperature and 0.1 MPa. Ref. [1, pp. 11 and 187] 



TABLE VIII: Calculated bulk moduli and elastic constants of maghemite. Unit: GPa. 





PW 


LCAO 




FiM GGA+C/ 


FiM GGA+C/ 


Bo 


146.3 147.8 


134.4 145.9 


Cll 


264.3 ± 27.0 285.0 ± 20.8 


245.2 + 12.3 266.3 + 31.2 


Cl2 


122.5 + 16.8 120.0 + 12.6 


114.2 + 9.0 113.2 + 22.7 


Cl3 


124.4 + 17.9 120.1 + 14.9 


113.1 + 1.1 114.2 + 15.4 


C33 


265.7 + 10.7 284.1 + 9.4 


246.0 + 4.3 266.7 + 12.7 


C44 


103.7 + 0.2 106.0 + 0.1 


90.9 + 2.7 94.3 + 3.2 


C55 


103.0 + 0.2 106.5 + 0.0 


92.4 + 1.8 95.8 + 6.5 



find that in all the calculations, cn ~ C33, cu ~ C13, and C44 ~ C55, which are conditions 
characteristic of the elasticity tensors of cubic crystals. This in indicative of the similarity 
between the tetragonal lattice of maghemite with its cubic counterpart. Although the long- 
range ordering of vacancies changes the symmetry of lattice, the elasticity tensor seems to 
be insulated from the change of symmetry. 
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E. Goethite (a-FeOOH) 



Goethite is the most thermodynamically stable iron oxyhydroxide, and has orthorhombic 
structure (space group Pnma, No. 62). 9l|] The lattice parameters have been measured 
by synchrotron powder diffraction at temperatures between 298 K and 429 K, j9l|, and at 



pressures up to 9 GPa, [92j. Gleason et al. |93| studied the equation of state of goethite 
under pressures 0-250 GPa, and found the equilibrium volume is 138.75 ± 0.02 A^, bulk 
modulus is 140.3±3.7GPa, and pressure derivative is 4.6±0.4. Goethite is antiferromagnetic 
in its ground state, with edge-sharing octahedron within a double-chain have antiparallel 
spin-moments, and corner-sharing octahedron in two double-chains have antiparallel spin- 
moments (see Figure E]). In addition to this antiferromagnetic state, we included another 
two antiferromagnetic states; one has same spin in a double-chain (noted as AFM'), the 
other one is similar to AFM except the corner-sharing octahedron have parallel spin (noted 
as AFM"). We have also calculated a ferrimagnetic, a ferromagnetic, and a non- magnetic 
state. 

The energetic order of antiferromagnetic states are the same in both PW and LCAO 
calculations. The energy difference between AFM and AFM' is only about 3 kJmol"^, 
which is, however, near the limits of the computation accuracy. The small energy difference 
between AFM and AFM' is reproduced in both the GGA+f/ and GGA calculations, using 
both the PW and LCAO basis sets. Despite this reproducibility across different basis sets, 
further calculations with high accuracy are required to distinguish the energetic order of 
the two antiferromagnetic states. The energies of AFM and AFM' are lower than AFM" 
by about 30 kJmol"^ in both PW and LCAO, indicating that corner- sharing octahedron of 
antiparallel spins (as in AFM and AFM') are energetically more stable than that of parallel 
spins (as in AFM"). In this study, we assume the AFM state is more energetically stable 
than AFM', and perform calculations of elastic properties based on the AFM magnetization 
state with and without on-site interaction. 

As we see from the calculations results in previous parts (hematite and maghemite), 
the lattice parameters from GGA+^7 are usually larger than that of GGA. This trend is 
violated in the calculation to lattice parameter c. In the PW calculations, the GGA+U 
result is smaller than GGA, while in the LCAO results, the calculated c are almost the same 
(Table HXl) . This feature is also seen in the calculation results of 6-axis of lepidocrocite (see 
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(b) AFM' 



^ * / ^ 

(c) AFM" 

FIG. 5: Antiferromagnetic configurations of goetliite. Cyan and blue octaliedron are for Fe 
atoms with antiparallel spin moments. Red balls are oxygen, white balls are hydrogen, 
sticks are hydroxyl bonds. Viewed along [010] direction. 

Section [III F[) . Since the hydrogen bonds are almost along c-axis in goethite (and 6-axis in 
lepidocrocite), the smaller values of lattice parameter c in goethite (and b in lepidocrocite) 
indicates strengthening of hydrogen bonds in the GGA+f/ calculations compared with the 
GGA calculations. The physical origin of this observation is not clear yet. Possible expla- 
nations may be from the redistribution of the charge density caused by the onsite Coulomb 
interactions. Although Fe atoms are not part of the hydrogen bonds (H-O..H), they have an 
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influence on the strength of hydrogen bonds by modifying the electron density in Fe-0 bonds 
which (in turn) change the electron density around the oxygen atoms, which are acceptors 
of the hydrogen bonds. The onsite Coulomb repulsion among Fe 3d electrons decreases the 
charge density in Fe-0 bonds, increasing the electron density around the oxygen atoms and 
strengthening the hydrogen bonds. The changes in electron density will be illustrated in 
more details in a separate paper. 

The ferrimagnetic state has the same average spin-polarization moment in both calcu- 
lations. Like the results of hematite, the spin moment of the ferromagnetic state is quite 
different: low spin in PW, and high spin in LCAO. With the exception of this difference in 
spin moments of ferromagnetic state, the calculation results from the PW and LCAO basis 
sets are consistent with each other. 

The calculated bulk moduli and elastic constants of the AFM state are listed in Table |Xl 
We found GGA+U calculations produce appreciably larger values of bulk moduli and most 
tensor components than the GGA in both PW and LCAO basis sets. The reason for the 
strengthening effect of GGA+U is not clear. It may be related to the hydrogen bonds 
which are sensitive to the distribution of electron density, but further work will be needed 
to understand this definitively. 



F. Lepidocrocite (7-FeOOH) 



Lepidocrocite has an orthorhombic structure (space group Cmc2i, No. 36 9J]), which 
consists of double chains of Fe(0,0H)6 octahedron which are aligned perpendicular to b- 
axis. The double chains form sheets, held together mainly by hydrogen bonds, which are 
weaker than covalent or metallic bonds, and may be longer than normal chemical bonds. 
Depending on the position of hydrogen atoms, the crystal structure of lepidocrocite can 



either be in the Cmcm space group (No. 63^ 
of two oxygen atoms in a hydrogen bond 



where the hydrogen atom reside at the middle 
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96| or in the Cmc2i space group (No. 36) 



where the hydrogen atom is closer to one of the two oxygen atoms. 9J] The difference is that 
Cmc2i is non-centrosymmetric, but is indistinguishable from the centrosymmetric Cmcm in 
XRD or neutron diffraction. The bond distances in the H-bonds in the Cmcm space group 
are extraordinarily large, thus the positions of hydrogen atoms may be averaged positions 
in neutron diffraction 96|, and the true symmetry may be Cmc2i (which has normal bond 
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TABLE IX: Calculated thermodynamic properties of goethite. See Table I for annotations. 





a (A) 


b(A) 


c 


;a) 


M(/iB) 


Ei{kJ mol-i) 


AFM 


10.018 (+0.6) 


3.017 


(-0.1) 


4.661 


(+1.2) 


0.00 


-453.1 


AFM' 


10.025 (+0.7) 


3.015 


(-0.2) 


4.650 


(+0.9) 


0.00 


-450.4 


AFM" 


10.045 (+0.9) 


3.055 


(+1-1) 


4.655 


(+1.0) 


0.00 


-419.6 


PW FiM 


10.039 (+0.8) 


3.042 


(+0.7) 


4.664 


(+1.2) 


2.50 


-432.9 


FoM 


10.103 (+1.5) 


2.842 


(-5.9) 


4.560 


(-1.0) 


2.63 


-417.6 


NM 


9.529 (-4.3) 


2.920 


(-3.3) 


4.366 


(-5.2) 


- 


-415.6 


GGA+U 10.040 (+0.8) 


3.045 


(+0.8) 


4.628 


(+0.4) 


0.00 


- 


AFM 


10.148 (+1.9) 


3.060 


(+1.3) 


4.654 


(+1.0) 


0.00 


-492.9 


AFM' 


10.149 (+1.9) 


3.061 


(+1.3) 


4.656 


(+1.0) 


0.00 


-490.0 


AFM" 


10.171 (+2.2) 


3.092 


(+2.3) 


4.683 


(+1.6) 


0.00 


-463.7 


LCAO FiM 


10.177 (+2.2) 


3.079 


(+1.9) 


4.673 


(+1.4) 


2.50 


-475.2 


FoM 


10.200 (+2.4) 


3.100 


(+2.6) 


4.695 


(+1.9) 


4.99 


-456.6 


NM 


9.642 (-3.2) 


2.939 


(-2.7) 


4.353 


(-5.5) 




-428.4 


GGA+C/ 


10.207(+2.5) 


3.105 


;+2.8) 


4.663 


(+1.2) 


0.00 




Expt.^'^) 


9.956 


3.021 


4.608 


0.00 


-559.3, -562.9 



(a) Measured at room temperature and 0.1 MPa. Ref. [1, pp. 11 and 187]. 



distances). We adopted the proposal in [9J] as the starting structure for our calculations. 

Each primitive cell contains 2 iron atoms, whose spin moment may align in parallel (fer- 
romagnetic) or antiparallel (antiferromagnetic) configurations. More magnetization states 
may be included if the magnetization state is stated in a conventional cell which contains 
4 iron atoms. Lepidocrocite is antiferromagnetic with aritiparallel spins in the same double 
layer, and antiparallel spins linked by hydrogen bonds. 94] This antiferromagnetic state is 
noted as AFM in this paper (Figure [HD. Another two antiferromagnetic states, noted as 
AFM' and AFM" are also included for comparison in addition to one ferrimagnetic, ferro- 
magnetic, and non-magnetic states. AFM' is similar to AFM, except the octahedron linked 
by hydrogen bonds have parallel spin moments; AFM" has parallel spin in a double layer, 
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TABLE X: Calculated bulk moduli and elastic constants of goethite. Unit: GPa. 





PW 


LCAO 




AFM 


GGA+U 


AFM 


GGA+U 


Bo 


93.1 


114.1 


98.6 


109.4 


Cll 


235.4 ±3.8 


298.4 ±4.4 


231.7 ±8.4 


252.2 ±5.3 


Cl2 


89.0 ± 1.2 


106.2 ±4.8 


86.6 ± 5.4 


89.2 ±0.4 


Cl3 


112.2 ±5.2 


117.4 ±3.6 


111.4 ±8.5 


117.7 ±6.0 


C22 


263.8 ± 4.3 


347.9 ± 3.6 


234.0 ± 10.0 


271.7 ±1.4 


C23 


96.3 ±6.1 


105.7 ± 2.0 


87.3 ± 11.1 


99.0 ± 1.7 


C33 


406.7 ±6.5 


414.8 ±3.9 


363.3 ± 12.0 


369.2 ±5.1 


C44 


78.9 ±0.1 


105.0 ±0.3 


58.7 ±0.2 


91.0 ±0.4 


C55 


122.6 ±0.7 


131.4 ±0.0 


106.9 ±2.8 


122.4 ±2.1 


C66 


65.9 ±0.0 


97.3 ±0.6 


60.1 ±2.1 


72.3 ± 4.8 



and antiparallel spin in neighboring double layers. 

The calculated thermodynamic properties of lepidocrocite are listed in Table Kit where 
we can see that calculations performed using the LCAO basis set produce larger error with 
respect to available experimental values than those obtained from the PW calculations. 
The largest error in our LCAO calculations is the overestimation to lattice parameter c 
by about 7.6%. The magnetization state of AFM' also deviated from antiferromagnetic, 
and converged to ferrimagnetic state during geometry optimization in LCAO calculations. 
The non-magnetic state has larger errors in both calculations than other states, which is 
consistent with the results obtained for the other iron oxides and oxyhydroxides, as described 
in previous sections. 

In this case, the energetic order predicted by PW and LCAO calculations compare very 
well among all the magnetization states, with the exception of the AFM' state (which devi- 
ates from the initial antiferromagnetic state) in LCAO calculations. The energy difference 
between the AFM and AFM" states are almost the same, about 10 kJmol"^, using both 
PW and LCAO. This indicates that the spin moments in a double-layer is unlikely to be 
parallel, as in the energetically unstable AFM" . The AFM" and FoM states have almost the 
same formation energies in both calculations, indicating there is weak correlation between 
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(a) AFM (b) AFM' (c) AFM" 

FIG. 6: Three antiferromagnetic configurations of lepidocrocite. Cyan and blue octaliedron 
are for Fe atoms with antiparallel spin moments. Red balls are oxygen, white balls are 
hydrogen, sticks are hydroxyl bonds. Viewed along [010] direction. 

iron atoms connected by hydrogen bonds. This can also be seen from a comparison of the 
AFM and AFM' states, between which the difference is solely due to the alignment of spin 
moments of iron atoms linked by hydrogen bonds. The formations energies of lepidocrocite 
(Table |XI|) are higher than those of goethite (Table HXl) . which is in good agreement with the 
established thermodynamic stability of the two iron oxyhydroxide phases (where goethite is 
known to be more stable than lepidocrocite). 

At this point we would like to highlight that a correct description of the loose, layered 
structure of lepidocrocite is much more difficult to obtain in our computations (using both 
PW and LCAO basis sets) than other types of oxides and oxyhydroxides. Geometry op- 
timizations often become trapped in an incorrect structure, as shown in Fig. [71 In the 
incorrectly optimized structure, the iron atoms are trigonal-bipyramid coordinated instead 
of octahedron, while the oxygen atoms that do not form hydroxyl bonds are bonded to only 
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TABLE XI: Calculated thermodynamic properties of lepidocrocite. See Table I for 

annotations. 



a (A) b(A) c(A) M(/iB) (kJ mol-i) 



AFM 


3.038 


(-1.4) 


12.624 


(+1.0) 


3.896 


(+0.7) 


0.00 


-425.8 


AFM' 


3.046 


(-1.1) 


12.604 


(+0.8) 


3.908 


(+1.0) 


0.00 


-425.2 


AFM" 


3.080 


(+0.0) 


12.204 


(-2.4) 


3.866 


(-0.1) 


0.00 


-417.4 


PW FiM 


3.066 


(-0.4) 


12.390 


(-0.9) 


3.898 


(+0.7) 


2.03 


-421.9 


FoM 


3.086 


(+0.2) 


12.233 


(-2.1) 


3.862 


(-0.2) 


4.17 


-416.9 


NM 


2.900 


(-5.9) 


11.846 


(-5.2) 


3.795 


(-1.9) 


- 


-398.2 


GGA+C/ 3.074 


(-0.2) 


12.546 


(+0.4) 


3.935 


(+1.7) 


0.00 


- 


AFM 


3.061 


(-0.6) 


12.417 


(-0.7) 


4.165 


(+7.6) 


0.00 


-463.5 


AFM' 


3.107 


(+0.9) 


12.533 


(+0.3) 


4.029 


(+4.1) 


2.47 


-457.8 


AFM" 


3.125 


(+1.5) 


12.456 


(-0.4) 


4.021 


(+3.9) 


0.00 


-452.7 


LCAO FiM 


3.107 


(+0.9) 


12.532 


(+0.3) 


4.029 


(+4.1) 


2.48 


-457.8 


FoM 


3.128 


(+1.5) 


12.472 


(-0.2) 


4.021 


(+3.9) 


4.94 


-452.0 


NM 


2.912 


(-5.5) 


11.745 


(-6.0) 


3.853 


(-0.4) 




-406.5 


GGA+C/ 3.091 


(+0.4) 


12.487 


(-0.1) 


4.011 


(+3.6) 


0.00 




Expt.^'') 


3.08 


12.50 


3.87 


0.00 


-554.6 



(a) Measured at room temperature and 0.1 MPa. Ref. [1, pp. 11 and 187] 



3 iron atoms instead of 4. The incorrect structure may have an abnormally small lattice 
parameter a (up to about 20% below experimental value), large b (up to about 25% above 
experimental value), or large c (up to about 40% above experimental value). These incorrect 
structures occurred when the geometry optimizations began using the structure models that 
have hydrogen atoms equidistant between oxygen atoms, and may also occur with certain 
computational settings. We show the equation of state (EOS) calculated using GGA and 
GGA+U to demonstrate the sensitivity of geometry optimization by GGA on the starting 
structure (Fig. [H]). The incorrect structure is accompanied by the steep energy decrease 
when the cell volume is slightly larger (3%; the correct structure can retain up to 1% vol- 
ume increase) than the equilibrium volume. In contrast, the GGA+f/ is robust in geometry 
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(a) Correct (b) Incorrect 

FIG. 7: Structural abnormality in geometry optimization. Gray blue balls are Fe, red are 
O, and white are hydrogen. Viewed along [010] direction, (color online) 

optimizations with varying volume in this case. We carefully examine the final structures 
after geometry optimization, and rigorously tested the computational settings and initial 
structures to ensure the double layer structure of lepidocrocite. 

In the case of lepidocrocite, the calculation results of elastic properties are more diverged 
than for other iron oxides and oxyhydroxides in this paper. As shown in Table IXIIt the 
difference between GGA and GGA+U can be more than 50% (cis, C22 and cqq in PW, C55 
and in LCAO), and the agreements between PW and LCAO are acceptable only for 
components C22, C23, and C33. This is partly due to the delicacy of lepidocrocite. 

G. Magnetite (Fe304) 

Magnetite has a cubic inverse spinel structure (space group FdSm, No. 227) at ther- 
modynamic standard state (room temperature, ambient pressure). Its chemical formula, 
Fe304, is often written as Fe'^"'"[Fe^"'", Fe^^]04 to show that tetrahedral sites (A) are oc- 
cupied by trivalent Fe ions, and octahedral sites (B) are occupied by equal trivalent and 
divalent ions. The spin moments of A- and B-sites align antiparallel, resulting in a ferri- 
magnetic state. Magnetite undergoes the Verwey phase transition at about 125 K, below 
which the electronic resistivity increases 2 orders of magnitude. jsTj] This phenomenon was 
explained by charge-ordering model in which electron hopping among Fe ions are frozen 
below Verwey transition temperature and aligned in an ordered pattern. |971] However, after 
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FIG. 8: Equation of state of lepidocrocite in GGA and GGA+?7 calculations, both using 
the PW basis set. The two vertical dashed lines mark the equilibrium volumes optimized 
by GGA and GGA+^7. The total energies are shifted to align the minimum energies of 

GGA and GGA+U. 



6 decades of study researchers found the phenomenon is far more complicated than was 
previously thought. [98] Among various changes (electronic resistivity, band structure, heat 
capacity) accompanied by the Verwey transition, the structure slightly distorted from the 
room-temperature cubic structure. At low temperatures, the structure of r nagnetit e was 



proposed to be orthorhombic from nuc 
oclinic from x-ray diffraction 



14, 



ear magnetic resonance spectroscopy 



99 



lOli mon- 



102|, neutron diffraction, [15|, elect ron d iffraction, [l^ 
and x-ray resonant scattering; [103l] or even lower symmetry of triclinic. 104| In the present 
study, we have restricted our calculation to the room-temperature cubic structure because 
the calculated thermodynamic properties at ground state can be extrapolated to room tem- 
perature without discontinuity by the phase transition. 

In the case of magnetite we have tested 3 magnetization states: ferrimagnetic (FiM), 
ferromagnetic (FoM), and non- magnetic (NM). Both PW and LCAO basis sets predict that 
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TABLE XII: Calculated bulk moduli and elastic constants of lepidocrocite. Unit: GPa. 





PW 


LCAO 




AFM 


GGA+C/ 


AFM 


GGA+U 


Bo 


75.6 


74.8 


X 


81.8 


Cll 


219.0 ±22.8 


246.9 ± 9.8 


242.4 ±0.9 


264.5 ± 12.8 


Cl2 


77.3 ± 14.5 


84.8 ±6.8 


88.0 ±0.2 


95.1 ± 13.2 


Cl3 


31.9 ±26.8 


80.6 ± 10.1 


79.0 ± 1.2 


72.2 ±3.8 


C22 


221.0 ± 11.6 


272.2 ± 18.7 


214.7 ± 1.7 


246.3 ± 10.4 


C23 


123.3 ±8.5 


137.2 ±22.9 


127.3 ±0.9 


124.1 ±4.9 


C33 


305.4 ±21.9 


347.7 ± 16.8 


303.2 ± 1.6 


327.4 ± 4.5 


C44 


121.6 ±0.0 


131.4 ±0.0 


97.6 ± 1.2 


119.0 ±0.8 


C55 


63.4 ±0.0 


64.0 ±0.1 


49.2 ±0.8 


63.6 ±1.9 


C66 


44.8 ± 0.0 


88.9 ±0.1 


73.3 ±0.6 


93.6 ±1.5 



the FiM state has the lowest formation enthalpy among all these magnetization states (see 
Table IXIIip . The lattice constants of the FiM state also provide better agreement with 
the experimental measurements. Using both basis sets, the calculation results differ from 
experimental values if we ignore spin polarization. The error in lattice constants using the 
LCAO approach is slightly larger than that using the PW approach. The calculated spin 
moments agree well with each other in both PW- and LCAO-based methods. 

We calculated elastic properties (bulk moduli and elasticity tensor) of FiM magnetite, 
as shown in Table IXIVI GGA calculations using PW and LCAO both predict that cubic 
magnetite is elastically stable. In the case of GGA+U calculations, we were unable to 
fit the strain energies with strains to calculate elastic constants, because the equilibrium 
cubic structure has higher energy than strained states. As shown in Figure [H except the 
isotropic deformation, the other two deformations have even lower energy than the zero- 
strain, "equilibrium" structure in GGA+U calculations. Increasing fc-point sampling density 
does not solve this problem. This indicates the cubic magnetite is elastically unstable, which 
agrees well with experimental observations that low-temperature (below Verwey transition 
temperature) structure is monoclinic, but not cubic. 
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TABLE XIII: Calculated thermodynamic properties of magnetite. See Table I for 

annotations. 







a (A) 


M(/xb) 


-E;f(kJ mol-1) 




FiM 


8.392 (-0.0) 


1.33 


-871.7 


PW 


FoM 


8.528 (+1.6) 


4.57 


-747.1 




NM 


8.049 (-4.1) 




-666.4 




GGA+C/ 


8.481 (+1.0) 


1.33 






FiM 


8.504 (+1.3) 


1.33 


-930.9 


LCAO 


FoM 


8.645 (+3.0) 


4.67 


-828.7 




NM 


8.111 (-3.3) 




-650.2 




GGA+C/ 


8.653 (+3.1) 


1.33 





Expt.(") 8.396 "-1120 



(a) Measured at room temperature and 0.1 MPa. Ref. |1, pp. 11 and 187] 

TABLE XIV: Calculated bulk moduli and elastic constants of magnetite. Unit: GPa. The 
calculation results with GGA+C/ are discussed in main text. 



PW 


LCAO 


FiM GGA+C/ 


FiM GGA+C/ 


B 


187.4 173.3 


165.3 168.3 


Cll 


275.4 + 40.9 


253.6 + 5.7 


Cl2 


155.2 + 60.3 


128.1 + 10.3 


C44 


97.5 + 13.0 


75.4 + 0.9 



H. Energetic order 

As mentioned above, one of the computational challenges in modeling different iron oxides 
and oxyhydroxides is the small energy differences among different solid phases. Since the 
typical accuracy of DFT calculations is about several kJmol"^, which is comparable to the 
energy differences between competing phases, calculations with different settings may lead 
to very different energetic order, but may have very little physical meaning due to numerical 
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inconsistencies. Systematic and consistent calculations of all the phases are highly desirable 
to make comparisons among the different phases, as well as case-by-case comparisons with 
experiments. The calculations in this study have enabled us to assess the PW and LCAO 
basis sets, but also to make such comparisons for the first time. 

To begin with, we have calculated the formation enthalpies of the iron oxides and oxy- 



hydroxides wit 
mental values 



1 re spect to hematite and water or oxygen (Figure [T0|) . Available experi- 



105| are also shown in figure for comparison. Among the 5 iron oxides and 



oxyhydroxides, magnetite is typically excluded from the experiments because all the other 
four compounds (goethite, lepidocrocite, maghemite, and hematite) can be rewritten as 
hematite+xH20 + AHf [x = for maghemite, or 0.5 for goethite and lepidocrocite). In 
our work, we are able to include magnetite and compare the energetic order with respect to 
hematite and oxygen. In this way we can plot the energetic order of all the 5 iron oxides 
and oxyhydroxides together, keeping hematite plus balancing gases (H2O and O2 in their 
standard state) as the reference. 

At this point, it is prudent to point out that our DFT calculations correspond to the 
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ground state, while experiments were conducted at thermodynamic standard state. The in- 
fluence of temperature and pressure on formation energies of solids are much less than that 
of gases. We therefore include the connection energy, whi ch d efines the difference between 



energies at ground state and standard state, for the gases. |106[ | By this ab-initio thermody- 
namics scheme, we can extend the calculation results at ground state to finite temperature 
and pressures. In this case, we use the connection energies for the gases which have been 



calculated to solve tribochemistry problems. 



lOTj For solids, the available thermochemistry 



data enable one to calculate connection energy by integrating from K ground state to the 
thermodynamic standard state as 

Afx'iZ) = r C,dT - T r ^dT, (2) 

where A/i°(Tr) is the connection energy at standard pressure = 1 atm) and room temper- 
ature Tr = 298.15 K, Cp is the molar heat capacity. The heat capacity, enthalpy difference 
between room temperature a nd K, and entropy can be looked up in thermochemistry ta- 



bles, eg NIST-JANAF table. 108| Calculation of the connection energies of gases are usually 



done through reaction equilibrium with solids. For example, we used the reaction 

MgO + H2O ^ Mg(0H)2 (3) 
to calculate the connection energy of g as-phase water because the thermochemical data of 



MgO and Mg(0H)2 are available. 107| Once the connection energies for room temperature 
is known, the chemical potentials at other temperatures (and pressures for gases) can be 
readily calculated using thermodynamics as long as the heat capacity data are available for 
the temperature range. 

From Fig. [10] we see the enthalpy difference between maghemite and hematite is appar- 
ently underestimated (by about 5 kJ mol~^) when using GGA with PW basis set, while all 
other settings reproduce this energy difference well. For magnetite, the enthalpy calculated 
using GGA+[/ with PW is larger than others, but the energetic order is consistent in all 
the calculations. The energy difference between goethite and lepidocrocite is larger in GGA 
than that in GGA+[/, in both PW and LCAO calculations. It is worth noting that the 
corrections using connection energy is ineffective to change the relative energetic orders of 
hematite and maghemite, or goethite and lepidocrocite, because they have same chemical 
compositions. However, with the corrections of connection energies, the relative energetic 
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FIG. 10: Calculated formation enthalpies of iron oxides and oxyliydroxides with respect to 
hematite and water or oxygen (The formation enthalpies of goethite and lepidocrocite are 
with respect to hematite and water; magnetite is with respect to hematite and oxygen). 

Expt stands for experimental values. 



order between compounds with different composition may change, as we see the sub-figures 
in Fig. [TOl With the corrections, the difference between the state for calculations and exper- 
iments are approximately eliminated, enabling us to make fair comparisons. The calculated 
energetic order of lepidocrocite, hematite, and maghemite is very different in the 4 data sets, 
and we may conclude that GGA+f/ with PW implementation best matches experiment. 

The consistent computational settings across different iron oxides and oxyhydroxides offer 
us a number of significant advantages, one of which is that we are in a position to construct 
phase diagrams. For this purpose we have chosen to use the calculation results from GGA+U 
with the PW implementation, and compute the the free energy of formation of a compound 
FeO^Hj^ as: 

AG = Ai/ - I (^Ao, (T) + i?r In ^) - I (ah, (T) + RTln^^ (4) 

where AH is the formation energy at ground state, Ao2{T) and AH2(r) are connection 
energies at a certain temperature for O2 and H2, respectively, R is the gas constant. One 
can write the formation energies with respect to H2O and O2 by analogy. The connection 
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energies were calculated using thermochemistry data in a previous study, lOTj and AH 
are from the ground-state calculations in the present study. In a phase diagram, the phase 
boundaries determined from equation H] are straight lines in a phase diagram. 

Using Equation H] we have constructed two phase diagrams, both corresponding to room 
temperature (see Fig. [TT]) . The metastable phases of lepidocrocite and maghemite are not 
shown, as these are equilibrium phase diagrams. The two sub-figures refer to the same 
systems with respect to the chemical potentials of (a) H2 and O2 and (b) H2O and O2, 
respectively. One notices the extremely low partial pressure of oxygen required for the for- 
mation of magnetite instead of hematite. This means that, at room temperature, magnetite 
should form under oxygen-poor conditions; otherwise the more stable hematite phase should 
prevail in exogenous environments. This is compatible with the fact that most magneto- 
tactic bacteria that produce magnetite are either anaerobic or microaerobic. |1, p. 481-489] 
Magnetite is able to form from hematite at low temperatures in the presence of hydrazine, l|, 
p. 405-406] which removes dissolved oxygen in the solutions. 

The phase boundary between hematite and goethite has the same slope of water formation 
in Fig. Illaj therefore, in a phase diagram of Fe-H20-02, it is independent of chemical 
potential of H2O, as shown in Fig. Illbi The phase diagram shows that the free energy 
of goethite is lower than hematite at standard state, and this agrees with the calorimetry 



measurements. pi)5| In a wet environment, these phase diagrams predict that the formation 
of goethite will be more thermodynamically favorable than hematite; while dehydration (dry 
conditions) will cause goethite to transform into hematite given a suitable driving force. 



IV. COMPARISON OF EFFICIENCY 



36, 
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481 ] The atomic orbitals 



One of the advantages of LCAO basis sets is efficiency. [32 , 
used to expand the wave functions are very economic (in terms of number of orbitals per 
atom to achieve an accuracy) compared with the PW basis set. In the DZP (double-^ with 
one polarization orbital) scheme which is used in the present study, each Fe atom needs 
17 orbitals for the valence electrons, O needs 13, and H needs 3 orbitals. For a 4 x 4 x 4 
A;-point grid for hematite, the number of atomic orbitals is 8832. In order to achieve similar 
convergence in energy calculations, the PW basis set requires about 120000 PW's, which is 



about 15 times as that of LCAO basis set. The advantage of less orbitals will be even more 
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FIG. 11: Phase diagram of iron oxides and oxyhydroxides with partial pressures of gases at 
T = 298.15 K. Lepidocrocite is metastable with respect to goethite, and maghemite is 

metastable to hematite; they are excluded in this equilibrium phase diagram. Dash-dotted 
lines indicate standard pressure (1 atm); the dashed line in (a) indicates formation of 
water vapor from H2 and O2 with no energy gain or loss. Note the different scales of 
partial pressures of H2in (a) and H2O in (b). 



apparent if the computation cell has vacuum space, such as in surface calculations, since the 
LCAO's are centered at ions, and vacuum requires no additional orbitals. In contrast, the 
PW's are delocalized, and even vacuum space has similar number- density of PW's. 

In addition to this, the localized nature of LCAO's enables one to implement the order-N 
algorithms, which critically rely on localization of wave functions. In integrating over bands, 
the Fermi level needs to reside in the band gap, which should be large enough to cover the 
varying chemical potential. This is not true for metals and semiconductors with narrow 
band gaps, which include most iron oxides and oxyhydroxides. Therefore, studies on these 
systems are not able to benefit from the order-N algorithms. 

However, fewer numbers of orbitals should still translate into efficiency (of computation 
time and memory usage), even without order-N algorithms. We find this is true for memory 
usage, but it is not always true for computation time. As shown in Fig. [121 the PW basis 
set uses more memory than LCAO for calculations of all the iron oxides and oxyhydroxides 
included in our study (see main text). It is worth noting that the memory requirements 
also depends on parallelization, and the numbers are extracted from calculations using 8 
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FIG. 12: Memory requirement of PW and LCAO in geometry optimizations to the iron 

oxides. 
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FIG. 13: CPU time consumption in geometry optimizations to hematite. 



CPUs for all the iron oxides except maghemite, which uses 32 CPUs. The CPU time usage 
for geometry optimizations for different magnetization states of hematite (Fig. [T3|) shows 
the PW basis set may exceed LCAO in some geometry optimizations, even though the PW 
basis set uses much more orbitals. Other factors may affect the computation time, such as 
minimization path, so we have taken care to always start from the same structures, and use 
the same method (CG) and force convergence (0.005 eV/A) in moving atoms in order to 
minimise this effect. The PW and LCAO basis sets also differ in their use of symmetry (as 
described above), which leads to differences in the force calculations. In general we find that 
the difference in computation time is not as large as that in number of orbitals. By utilizing 
symmetrization. As an aside, we also compared the numbers of self-consistency iterations 
to reach the geometry optimization criteria. In most cases, the LCAO basis set needs more 
MD steps than the PW basis set to reach the convergence criteria. 
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At this point it is also worth pointing out that one of the problems with LCAO's is 
systematic convergence. Simply increasing radii of the atomic orbitals does not always lead 
to better convergence, and tuning the parameters of the atomic orbitals requires considerably 
more effort than is needed for the PW basis set. While increasing the number of atomic 
orbitals can increase the accuracy, this comes at the cost of computation (in the PW basis 
set as well) . Tests of the size of atomic orbitals have shown high accuracy within the frame 
work of DFT can be achieved with multiple-^ and multiple polarization orbitals. DZP, 
which is used in the present study, is usually a reasonable compromise between accuracy 
and efficiency. 

V. CONCLUSIONS AND DISCUSSIONS 

In summary, by comparing the calculation results of PW and LCAO basis sets, with and 
without on-site interactions, as well as among different magnetization configurations, we 
presented solutions to the computational challenge in modeling iron oxides. Consistency is 
paramount, and this has been maintained in the comparisons as to energy functionals, con- 
vergence criteria of force, /c-point mesh, and starting structures for geometry optimizations. 
We have shown that both PW and LCAO basis sets can find the thermodynamically stable 
magnetization states, and reproduce lattice parameters well (except lepidocrocite by LCAO 
which overestimate c by about 7% in GGA and 4% in GGA+U). However, in most geom- 
etry optimizations, the LCAO basis set is more efficient in CPU time and memory usage 
than the PW basis set, but the accuracy is slightly reduced when comparing with PW basis 
set. Several factors contribute to the efficiency difference between the two implementations, 
including number of orbitals, molecular dynamics algorithms in moving ions, electron den- 
sity mixing, force calculation, A;-point density. Using these basis sets, we evaluated elastic 
stability of all the materials. We find that the PW and LCAO basis sets are comparable 
for most structures except lepidocrocite, and that the elasticity tensor of maghemite is close 
to that of a cubic crystal, though the true symmetry is tetragonal due to the long-range 
ordering of vacancies. While GGA predicts cubic magnetite is elastically stable, GGA+U 
calculations contradict the prediction. 

The crystal structure of lepidocrocite consists layers held by H-bonds. In computational 
modeling, functionals with general gradient approximation and hybrid functionals are able 
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to describe the relatively weak interactions of hydrogen bonds. While van der Waals inter- 
actions may also contribute significantly to inter-layer interactions, as they do in graphite, 
they are not included in the present study, as calculations of dispersive forces are either very 
expensive or relying on empirical parameters. In our calculations this delicate structure ex- 
hibited some structural abnormalities, which may be due to the omission of dispersive forces. 
This layered structure is not as delicate as that of graphite (in which the carbon layers are 
held by even weaker van der Waals interaction), but still impose an challenge to computa- 
tional modeling. Accurate energy functionals that include van der Waals interactions may 
describe better the crystal structure of lepidocrocite. 

Based on these results (which represent the first consistent set of ab initio predictions 
of the clastic, magnetic and thermodynamic properties), we also present the first phase 
diagram of 5 iron oxides and oxyhydroxides designed to predict the relative stability of these 
materials under different chemical conditions. Given that chemical conditions are typically 
characteristic of specific environments (both during and post- formation) , this phase diagram 
will be invaluable in understanding the environmental stability of these important materials, 
and anticipating transformations that may be invoked by moving from one environment to 
another, or by variations in climatic conditions. 
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